A Simple method to set up low eccentricity initial data for moving puncture 

simulations 
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We introduce two new eccentricity measures to analyze numerical simulations. Unlike earlier 
definitions these eccentricity measures do not involve any free parameters which makes them easy to 
use. We show how relatively inexpensive grid setups can be used to estimate the eccentricity during 
the early inspiral phase. Furthermore, we compare standard puncture data and post-Newtonian 
data in ADMTT gauge. We find that both use different coordinates. Thus low eccentricity initial 
momentum parameters for a certain separation measured in ADMTT coordinates are hard to use 
in puncture data, because it is not known how the separation in puncture coordinates is related to 
the separation in ADMTT coordinates. As a remedy we provide a simple approach which allows 
us to iterate the momentum parameters until our numerical simulations result in acceptably low 
eccentricities. 

PACS numbers: 04.25.dg, 04.25.Nx, 04.20.Ex, 04.30.Db, 



I. INTRODUCTION 



Currently several gravitational wave detectors such as 
LIGO [EEf, Virgo @, i] or GEO @] are already operat- 
ing, while several others are in the planning or construc- 
tion phase [6]. One of the most promising sources for 
these detectors are the inspirals and mergers of binary 
black holes. In order to make predictions about the fi- 
nal phase of such inspirals and mergers, fully non-linear 
numerical simulations of the Einstein Equations are re- 
quired. To numerically evolve the Einstein equations, at 
least two ingredients are necessary. First we need a spe- 
cific formulation of the evolution equations. And second, 
to start such simulations initial data are needed. As the 
first ingredient most groups nowadays use the BSSNOK 
formulation 0-Q of the evolution equations. This formu- 
lation is usually evolved using finite differencing methods, 
but for single black holes there have been some attempts 
to use spectral methods (Iol - [i2| . For binary black holes 
the BSSNOK system is usually used together with the 
moving puncture approach [l2| [l4j]. This approach so 
far works only with finite differencing methods since cer- 
tain evolved variables are not smooth inside the black 
holes (at the punctures). Almost all simulations using 
the BSSNOK formulation to date use standard puncture 
data [HI, [l6[ as initial data. These initial data are very 
flexible in that they contain free parameters for the posi- 
tion, momentum and spin of each black hole and thus one 
can setup practically any kind of orbit. Note, however, 
that the emission of gravitational waves tends to circu- 
larize the orbits [f7|, EH ■ Thus for realistic binary black 
hole systems that have been inspiraling already for a long 
time, we expect the two black holes to be in quasi-circular 
orbits around each other with a radius which shrinks on 
a timescale much larger than the orbital timescale. This 
means that the initial data should be such that the or- 
bit has no or at least very small eccentricity. For our 
purposes here we follow the NRAR (Numerical Relativ- 
ity - Analytical Relativity) collaboration [l9[ guidelines 



which consider eccentricities of order a few times 1 0~ 3 ac- 
ceptably small. There have been several previous works 
that have considered eccentricities for puncture initial 
data [2fj| - |22j . However, the most successful approach in 
terms of achieving low eccentricities was implemented for 
excision type initial data [231-25] . The method discussed 
in this work aims at lowering the eccentricity for the kind 
of puncture initial data that is routinely used with the 
moving puncture approach. 

Throughout we will use units where G = c = 1. The 
black hole masses are denoted by mi and m^. We also 
introduce the total mass M = mi+1112, the reduced mass 
/i = m\milM and v = fi/M. 

The paper is organized as follows. Sec. |TT] introduces 
and compares several eccentricity measures. In Sec. IIIII 
we describe grid setups that can be used in numerical 
simulations aimed at measuring the eccentricity. Sec. IIVI 
discusses a simple method to pick initial momentum pa- 
rameters. This is followed by Sec. IVl which describes how 
to iterate these parameters to arrive at a reduced eccen- 
tricity. We conclude with a discussion of our results in 
SeclVH 



II. 



DEFINING ECCENTRICITY FOR INSPIRAL 
ORBITS 



Real binary black hole orbits can never be circular. 
They always follow spirals. So when we are aiming for low 
eccentricity initial data, we really want data that result in 
trajectories that spiral in smoothly without oscillations in 
the black hole separation. Of course this issue is further 
complicated by the fact that trajectories are coordinate 
dependent. 

There are several earlier eccentricity definitions for in- 
spiral orbits in the literature [20i-25]. All of them define 
eccentricity as a deviation from an underlying smooth, 
secular trend in some _specific quantity that is associated 
with the orbits. In [20| the frequency of the dominant 
I = m = 2 mode of the gravitational waves emitted is fit- 
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ted to a fourth order monotonic polynomial, and the de- 
viation of the frequency from this fit is used to compute 
the eccentricity. This approach works for non-spinning 
binaries. Essentially the same method is also used in (2l| . 
but instead of the gravitational wave frequency [2lj uses 
the orbital frequency and also the coordinate separation 
to obtain two eccentricity measures. These same mea- 
sures are also used in [22[ . The approaches in [23|, [24[ fit 
a linear function plus a sine function to the coordinate 
separation and also the proper separation. The eccen- 
tricity can then be obtained from the amplitude of the 
fitted sine function. In [25j j the same fitting approach as 
in [21| is used, but the fitted quantities are coordinate 
separation, proper separation and also orbital frequency. 
All the approaches based in orbital parameters should in 
principle also work for systems with spin. Also note that 
all these eccentricity definitions are chosen such that they 
result in the correct value for Newtonian orbits. 

Below we will introduce two new eccentricity defini- 
tions and compare them to the earlier definition based 
on fitting the orbital frequency [2l], [25[ . 

The first eccentricity definition is based on the coordi- 
nate separation of the two black holes. It is given by 
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where the average separation, and the maximum and 
minimum deviation from a smoothed value r s are given 
by 

rt+T/2 

r av = / r(t')dt'/T (2) 

Jt-T/2 

Ar max (t) = max [r(t') - r s {t', t)} (3) 

t'e[t-T/2,t+T/2Y 

Ar min (t) = min [r(t') - r s (t', <)]. (4) 

t'e[t-T/2,t+T/2Y 

Here the period T is defined using Kepler's law 

T = 2vr(r 3 /M) 1/2 . (5) 

Notice that the actual orbital period may be slightly dif- 
ferent, but this estimate suffices to get an approximate 
eccentricity measure. The smoothed value r s (t',t) is ob- 
tained from 



r s (t',t) = r{t) 



r(t + T/2) - r(t - T/2) 
T 



(t'-t), (6) 



but different smoothings are possible (e.g., by performing 
a least-squares fit of e.g. a linear or quadratic function 
to r(t) in the interval [t + T/2, t - T/2]). Essentially the 
definition in Eq. ([1]) measures how much the coordinate 
separation oscillates over the time T. For Newtonian 
orbits it coincides with the usual eccentricity definition 
for elliptic orbits. For orbits whose radius shrinks linearly 
in time (without any oscillations) e r {t) is zero. 

Another similar eccentricity measure can be obtained 
using the gravitational wave signal of the inspiraling 



binary. The idea is to determine the separation in a 
more gauge invariant way from the amplitude of ^4 in- 
stead of using the gauge dependent coordinate separa- 
tion. In [26j j it is shown that for a non-precessing binary 
in the quadrupole approximation the amplitude of the 
I = m = 2 spin-weighted Spherical Harmonic mode is 
given by 



\C 2 



32v^~/5v(Mlu) 8/3 , 



(7) 



where lu is the orbital angular velocity. Using Kepler's 
law we can define a separation 



f22 



M i/3 w -2/ 3 = m[|C 22 |/(32 VV5^)]- 1/4 (8) 



which is directly related to the amplitude of IC22I of the 
I — Tfi — 2 mode of ^4. Replacing the coordinate sepa- 
ration in Eq. ([!]) by r 22 we define 



eaa(i) 



Ar 2 2,maa(Q - Ar 2 2,mmft) 
2r 2 2.av 



(9) 



which is an eccentricity definition that can be computed 
from \&4 alone. Note that this definition needs to be 
extended for the case of precessing orbits, since in that 
case |C 22 | will oscillate even for spherical orbits (i.e. or- 
bits with r = const). The extension could be achieved 
by instead using a |C 22 | that is computed in a coordinate 
system where the z-axis points along the instantaneous 
orbital angular momentum. 

We have also tested an eccentricity definition based on 
the coordinate angular velocity to. Here the eccentricity 
is defined by [H l2l| 



Uj(t) - UJfjtjt) 



(10) 



where u>(t) is simply the coordinate angular velocity and 
Wfit{t) is a polynomial fit of order 5 to cj(t) over a time 
interval corresponding to several complete orbits. The 
hope is that the fit will smooth out oscillations so that 
e u (t) oc uj(t) —ujfit(t) becomes a measure of how much uj 
oscillates. The actual eccentricity is the maximum mag- 
nitude of e^it). 

In order to compare these three eccentricity definitions 
we now present the eccentricities from an actual numer- 
ical simulation. As we can see in Fig. [1] all three eccen- 
tricity definitions agree well for 300M < t < 3000M. 
They yield a value of about 0.02. Notice that e r and e 22 
are direct eccentricity measures, while in the case of e u 
the eccentricity corresponds to the maximum value of the 
magnitude of e w . The eccentricity definition e r which is 
calculated from the separation has certain problems for 
t < T/2 ~ 135M. These occur because we have no data 
for t < 0, which is needed in the average over one com- 
plete period (centered around t) and also for the fitting 
in Eq. ([6]). These same problems also affect e 22 . In e 22 , 
however, they are exacerbated by the initial junk radi- 
ation that dominates |C 22 | until about 150M (see mid- 
dle panel). We see that e 22 does not completely settle 
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NOK formalism 0-0 in the variation known as the "mov- 
ing punctures" method [U, 13 ■ The particulars of our 
BSSNOK implementation can be found in [H H^. For 
completeness we note that lapse and shift evolve accord- 
ing to 



1500 2000 
t/M 

FIG. 1: This plot shows results from a numerical run (with 
parameters from row 6 of table [I] The upper panel shows 
the coordinate separation r. The middle panel depicts the 
magnitude of the / = m = 2 mode |C7 22 1 of ^4 extracted at a 
separation of 70Af from the center of mass. In the lower panel 
we plot the three different eccentricities obtained with the 
definitions given in Eqs. (JXJ) , © and (|10p . Note that e r and 
e r directly measure eccentricity, while for e u the eccentricity 
corresponds to the maximum values of |e„,(i)|. 



down until about 700M. Notice also that during stan- 
dard moving puncture evolutions the coordinates adjust 
quite rapidly initially. This means that any eccentricity 
definition that is based on the coordinate separation or 
the coordinate angular velocity will not be completely 
reliable during the first 100 M or so. Thus it comes as no 
surprise that the eccentricity definition e u has a different 
frequency and amplitude in the beginning. 

The curves for e r and e22 in Fig. Q] oscillate even at 
later times. This oscillation is due to the fact that the 
period T which we get from Kepler's law is not exactly 
equal to the actual orbital period. The magnitude of this 
oscillation can be used as an error estimate of our eccen- 
tricity measures. For a conservative eccentricity estimate 
we can use the maximum values of e r and e22- 

Let us point out that e r and e22 are easier to com- 
pute than e u . The latter depends on a polynomial fit 
to the measured orbital angular velocity. The problem 
is that this fit has to be done over a certain time inter- 
val, that must terminate well before the merger. Thus 
it requires a certain amount of fine tuning and human 
intervention. On the other hand e r and e22 can be com- 
puted at any time in a very simple way. Recall however, 
that e22 is more sensitive to the initial junk radiation. 
Thus for short runs used to probe the eccentricity of a 
configuration we usually just use e r . 



III. NUMERICAL EVOLUTIONS 

The numerical results discussed in this paper have been 
obtained with the BAM code As already men- 

tioned, the gravitational fields are evolved using the BSS- 



(fit - ?di)a 
(d t - p k d k )(3 l 



-2aK, 



{d t -P k d k )B l = (d t ~f3 k d k )f l ~ V B l . (11) 

The shift driver parameter is set 77 = 2/M in all our runs. 

The BAM code is based on a method of lines approach 
using sixth order finite differencing in space and explicit 
fourth order Runge-Kutta time stepping. The time step 
size is chosen such that the Courant factor is either 0.25 
or 0.5. For efficiency, Berger-Oliger type mesh refinement 
is used [30l | . The numerical domain is represented by a 
hierarchy of nested Cartesian boxes. The hierarchy con- 
sists of L + 1 levels of refinement, indexed by I = 0, . . . , L. 
A refinement level consists of one or two Cartesian boxes 
with a constant grid-spacing hi — h /2 l on level I. We 
have used here L = 10 to 11 for the number of refinement 
levels, with the levels through 5 each consisting of a sin- 
gle fixed box centered on the origin (the center of mass). 
On each of the finer levels 6 through L, we initially use 
two sets of moving boxes centered on each black hole. 
When the black holes get close enough that two of these 
boxes start touching, they are replaced by a single box. 
The position of each hole is tracked by integrating the 
shift vector. We use this same set up but with different 
resolutions depending on the purpose of each simulation. 

For an accurate simulation of the inspiral and merger 
of two non-spinning equal mass black holes we might use 
L = 10 with a resolution h%o = M/96 on the finest level 
using 144 points on the fixed levels and 72 points on the 
moving levels. The notation we use to describe this grid 
setup for this simulation is: 

[5 x 72,6 x lU][M/h w = 96, OB = 768M][C = 0.25] 

which indicates that we have 5 moving levels with 72 
points in each box and 6 fixed levels with 144 points 
each. The resolution is given by M/hio — 96 on the finest 
level, which results in an outer boundary at 768 M. The 
Courant factor here is chosen to be 0.25, which implies a 
time step of dt\o — 0.25/iio = Af/384 on the finest level. 
If the black holes have spins and / or unequal masses even 
more resolution is needed. For example for a mass ratio 
of 3 and a dimensionless spin magnitude of 0.6 on the 
larger hole we might use a setup described by: 

[6 x 72, 6 x 144] [M/h n = 192, OB = 768M] [C = 0.25] 

(13) 

this setup has twice the resolution on the finest level, so 
that the resolution in terms of the individual masses is 
now given by m\/hn — 48 for the smaller and rtiijhw = 
144 for the larger hole. Such runs are quite expensive. 
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FIG. 2: This plot shows the eccentricity e r for two different 
grid setups. Both start with the parameters from row 7 of 
table[I] The solid line shows e r for a high resolution grid, while 
the broken line shows the measured e r if we use a less accurate 
and computer intensive grid setup. Both yield similar orbital 
eccentricities. 



Until merger they take about one month on a Cray XT5 
supercomputer like NICS Kraken if we use 96 cores. 

However, if our objective is to simply measure the or- 
bital eccentricity of a binary (characterized by certain 
initial parameters), it is sufficient to evolve for only a 
few orbits. We have found that such an evolution does 
not need great accuracy. So if the goal is to simply de- 
termine the initial eccentricity we use the following setup 

[5 x 48, 6 x 96] [M/h w = 85.3, OB = 576M] [C = 0.5] 

With this grid setup it takes only two days to evolve up 
t = 800M if we use 48 cores on NICS Cray XT5 Kraken. 
In Fig. [5]we see that the eccentricity e r from this cheaper 
run agrees quite well with the e r from a more expensive 
run. As explained above no eccentricity estimate is ac- 
curate at the start of a run. We need to wait until about 
500M for e r to settle down to a regular oscillation. As 
mentioned before we use the maxima of e r as a conserva- 
tive eccentricity estimate. Thus the eccentricities of the 
cheap and expensive runs are 0.0013 and 0.0010. Hence 
we need to evolve until at least about 800M to get a reli- 
able estimate. From the oscillations in the curves we see 
that the errors in both eccentricity estimates are about 
0.0005. 



IV. CHOOSING PUNCTURE PARAMETERS 

In order to start our simulations we need initial data 
for binary black holes with arbitrary spins and masses at 
some given initial separation r. Since we will employ the 
moving punctures approach in our evolutions, we will use 
standard puncture initial data (l5j . Thus the 3- metric 
and extrinsic curvature are given by 



A=l 



pW 



\e-kimb A n A d y n A > 



8kiP k A nl 



A (gij 



n\n> A ) 



(16) 



Here p A and S\ are the momentum and spin parameters 
of black hole A, while r A and n l A denote the distance 
and normal vector measured from hole A. The conformal 
factor is 



jh — 1+ 1 H±L + H^i 
2n 2r 2 



(17) 



where nib 1 and mt 2 denote the black hole bare mass pa- 
rameters. The scalar u is computed by numerically solv- 
ing the Hamiltonian constraint. These initial data are 
very flexible since the free parameters for the position, 
momentum and spin of each black hole can be chosen 
freely. Thus one can setup practically any kind of orbit. 
Note, however, that our goal is to set up data for black 
holes that are in quasi-circular orbit. This means that 
we need to choose our momentum parameters such that 
the eccentricity is as small as possible. Since we start 
our evolutions in a frame where the center of mass is at 
rest, both black holes have momenta that are equal in 
magnitude but opposite in direction. Thus we have to 
choose only two parameters: the tangential and radial 
component of the momentum of one of the black holes. 

To complete the definition of the initial data, we also 
need to specify initial values for the lapse a and shift 
vector /?'. At time t = we use 

a = V" 2 , 

& = 0. (18) 

A. Finding the momentum parameters 

There have been various attempts to guess appropri- 
ate momentum parameters. Some have used the quasi- 
equilibrium approach [3lM33j , but most are based on 



34. 
21] 



9*3 



(15) 



post-Newtonian (PN) approximations (see e.g. [29j, 
|35|). The latter have been taken to their extreme in 
and [2a who integrate PN equations in the ADMTT 
gauge [36} in the hope that the thus obtained momen- 
tum parameters will lead to orbits with less eccentricity. 
However, as mentioned in [22j], it is not certain that non- 
eccentric PN parameters in ADMTT gauge should pro- 
duce non-eccentric orbits in full General Relativity. We 
want to point out here, that standard puncture data are 
inconsistent with PN theory beyond (u/c) 3 [46| (34l.l37r- 
[39j . even in ADMTT gauge! The reasons for this incon- 
sistency are as follows. First, the 3-metric of puncture 
data is always conformally flat (see Eq . (TTSlO . while the 
PN 3-metric contains deviations from conformal flatness 
at order (v/c) 4 [13, S3). This is true for both harmonic 
and ADMTT gauge. Second, the conformal factor in 
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ADMTT gauge is given by [13 

i> PN = l + ^ + ^, (19) 

where Ea is the energy of each particle used to model 
the black holes. If we compare Eqs. (|19p and p7|) we 
see that the conformal factor in ADMTT gauge is not 
identical to the conformal factor used in puncture initial 
data. One difference is that the ADMTT ippN contains 
the particle energies Ea while the puncture ip contains 
the bare masses m bA . These two only agree for infinite 
separation. Furthermore, the puncture ip contains the 
additional piece u that is obtained by numerically solv- 
ing the Hamiltonian constraint, while PN data violates 
the Hamiltonian constraint. One might like to argue that 
the PN data should agree with puncture data up to the 
higher order terms that are neglected in the PN approxi- 
mation. This is however not true close to the black holes, 
because PN theory is not valid in regions of very strong 
gravity. Thus near the black holes the ADMTT metric 
differs from the puncture metric by terms that are of low 
PN order. 

From the above explanation it is clear that the coor- 
dinate separation in ADMTT gauge does not have the 
same physical meaning as the coordinate separation in 
puncture initial data. Thus if by some method we arrive 
at the momentum parameters needed at a particular sep- 
aration r (in ADMTT coordinates), we should not simply 
use these same momentum parameters at the same punc- 
ture separation. 

So if (for some configurations) the more complicated 
approaches in 21 1 or [22| lead to less eccentric orbits than 
e.g. the simpler approaches in [2^, H3, [H[ this should 
be considered a coincidence, since it is not due to the 
inclusion of higher order terms or due to the integration 
of post-Newtonian equations of motion. In fact, we find 
that (as anticipated by |22]) for many configurations the 
approaches in [21 1 or 1221 do not lead to less eccentricity 
than the ones in [29|, [3J, T35| . 

For these reasons we take a different approach here. 
We will take a simple PN formula to obtain a reason- 
able guess for both the tangential and radial momentum 
components pt and p r . Then we numerically evolve the 
resulting initial data for a short time (using the efficient 
grid setup in Eq. (|14p) to see how eccentric they really 
are. Afterward we adjust p t to reduce the eccentricity. In 
this way we can obtain low eccentricity orbits for any con- 
figuration. In order to come up with a reasonable guess 
for pt a nd p r we use 2PN accurate expressions of Kid- 
der 4l| in harmonic gauge. Specifically, we freely choose 
the two masses mi, m.2, the six spin components of Si 
and S2 and a separation r. Next, choosing L jy in the in- 
direction, we use Eqs. (2.8) and (4.7) of [4lj to compute 
the total orbital angular momentum L, and Eq. (4.12) 
of [4l| to compute r. We rotate L, Si and S2 so that L 
points in the z-direction. Then we obtain the momentum 
in the xy-plane as 



p r = n\r\ 



(21) 



In all cases discussed in this paper we put the two punc- 
tures on the y-axis at y± = m 2 r/M and 1/2 = —m\rjM. 
The two initial black hole momenta are then pi = 
(-pt,—Pr,0) and p 2 = (pt,Pr,0). 



B. Determining the bare mass parameters 

The momenta pi, P2 and the spins Si, S2 for a co- 
ordinate separation r directly enter the Bowen-York ex- 
trinsic curvature of standard puncture data. Note, how- 
ever, that the bare mass parameters m bl and m& 2 which 
appear in the construction of standard puncture data 
are not equal to the individual black hole masses. As 
in [29|, [34| 5H, [43[ we obtain the bare masses from the 
condition that the ADM masses 



m 



ADM 



m bA (l + u A ) + 



m bl m b2 
27 ' 



(22) 



measured at the punctures [44| should be equal to mi 
and m 2 , where ua is the value of u at puncture A. As 
in [l6|, HH, Hi] we assume that the ADM masses measured 
at each puncture are a good approximation for the initial 
individual black hole masses. Numerically this condition 
is implemented as a root finder in the initial data solver 
that picks m bl and m b2 such that the ADM masses at 
the punctures are equal to mi and m.2 . 



V. REDUCING THE ECCENTRICITY 

Now that we know how to generate initial data on 
approximately circular orbits for arbitrary spins, masses 
and separations, it is time to present some numerical re- 
sults to demonstrate what eccentricities we get and how 
we can reduce them. For our purposes here, we consider 
the eccentricity small enough if it is 0.003 or less, which 
agrees with the NRAR target of an eccentricity of a few 
times 10 -3 . The first example we consider is an equal 
mass binary without spin. The momentum parameters 
are picked according to Eqs. (l20l) and (|2~Tj) . The values 
of all the initial parameters are given in the first line of 
table |U As we can see, the eccentricity is about 0.005 in 
this case. In order to test how it varies with p t we have 
also performed a run where we have increased pt by 



Ap t = p, 



, MY (M 
11.29 1— - 92.37 — 
r J \ r 



(23) 



Pt = |L|/r 



(20) 



The results of this increase are shown in the second line 
of table U and yield an eccentricity that is reduced by 
more than a factor of 3. The expression for Eq. (|2"3")l 
comes from the fitting of two of our older equal mass runs 
that resulted in low eccentricity. Thus non-zero spins or 
unequal masses are not taken into account by this fit. 
Therefore Eq. (f2"3"]l certainly does not present the opti- 
mum momentum correction. In fact we have found that 
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2/3 
2/3 
2/3 


0.4, 0, 

0.4,60,0 

0.4,60,0 

0.4,60,0 

0.3,0,0 


0.6, 180, 
0.6,60,0 
0.6,60,90 
0.6, 120,90 



6.9756 
7.6504 
7.6509 
7.7896 
7.7196 


3.19 
3.18 
3.19 
3.51 
3.32 


3.0 
2.0 
1.8 
2.1 
1.3 



TABLE I: Initial data parameters. The black holes have coordinate separation r. We give both bare masses mt 1 , mj 2 as well 
as physical masses mi, m,2. The punctures are located on the y-axis at yi = m,2r /M and y2 = —mir/M. The spins are given 
in terms of their magnitudes and the usual polar angles of spherical coordinates measured in degrees. The linear momenta are 
(-fPt, TPr, 0). The last column shows the resulting eccentricity. 



adding Ap t to pt does not always reduce the eccentricity 
since the PN estimate of Eq. (j2"0")) for p t is sometimes too 
large and sometimes too small for generic orbits. The 
expression in Eq. (|23l) simply gives us a rough estimate 
by how much we might have to raise or lower p t in order 
to reduce the eccentricity. 

So in order to really reduce the eccentricity we usually 
start one run on the coarse grid described by Eq. (TT4l . 
This run normally uses pt given by Eq. ([20]) . We then 
look at the coordinate separation r(t) for this run. Usu- 
ally one can tell whether the initial tangential momentum 
Pt was too large or too small. We then start a new run 
where we either increase or decrease p t by Ap t . This run 
then gives a new r(t) and e r . From the two resulting 
simulations one can then extrapolate to zero eccentricity 
to obtain a more refined tangential momentum parame- 
ter. This procedure is illustrated in rows 3,4, and 5 of 
table |U In row 3 we have used the p t = 8.5013 x 10 _2 M 
from Eq. ([20)1 and obtained an eccentricity of 0.004. In 
row 4 we have increased this pt by Apt which leads to a 
decrease of the eccentricity to 0.002. A further increase 
to p t = 8.5280 x 10~ 2 M (as in row 5) then yields an 
eccentricity of 0.0015. 

Row 6 of tableUshows the initial data that were used to 
produce Fig. [I] In order to produce an eccentricity which 
is clearly visible in the r(t) curve we have used a pt that 
is deliberately chosen much larger than what is predicted 
by Eq. (1201) . In row 7 we show the same configuration 
(at somewhat closer separation), but this time we choose 
Pt according to Eq. (|2"0"|) . we see that in this case the 
eccentricity is already small enough to satisfy e.g. the 
NRAR guidelines. 

In row 8 we have used our method to produce momen- 
tum parameters for a similar configuration, but this time 
both spins point in the negative z-direction. Again we 
can reach the NRAR target of an eccentricity of order 
0.001. In row 9 we also show the eccentricity resulting 
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FIG. 3: This plot shows r(t) for the two numerical runs with 
parameters from rows 8 (solid line) and 9 (broken line) of ta- 
ble U The parameters for the solid line were picked according 
to the method explained in this paper while the broken line 
is the result of setting the parameters according to the much 
more complicated method in [52] . As we can see the eccen- 
tricity difference of about a factor of 5 is quite noticeable in 
the r(t) curves. 



from the method introduced in [22| for the same mass 
ratio and spin configuration. We can see that it is about 
5 times larger, which supports our argument that inte- 
grating PN equations of motion does not necessarily lead 
to better results. The results from these two runs are 
also shown in Fig. [3] The solid line corresponds to row 
8 (i.e. our approach) and the broken line is generated 
using the approach in [22]. As we can see our method 
does give a noticeably less eccentric r(i) curve. Notice 
that the initial dip in the coordinate separation r{t), is 
due to the aforementioned initial coordinate adjustment. 
It does not mean the two holes really plunge toward each 
other initially. 
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The last five rows of tableQ]give a few more examples of 
parameters and resulting eccentricities. They show that 
we can obtain low eccentricities for generic mass ratios 
and spin orientations. 

Notice that our eccentricity reduction method is simi- 
lar in spirit to the method described in [23|, Hil , in that 
we also use short numerical runs to adjust certain pa- 
rameters. There are, however, important differences. 
In [23|, excision type initial data, not punctures are 
used. These data are constructed using an extension 
to the conformal thin sandwich formalism [45| . Empir- 
ically it turns out that the tangential momentum that 
is achieved with the conformal thin sandwich method is 
quite close to what is needed for low eccentricity data. 
The main reason for eccentricity is the absence of a ra- 
dial momentum in the original conformal thin sandwich 
method. In [23| a method is developed that adds an arbi- 
trary radial velocity parameter to the initial data. This 
radial velocity parameter together with the tangential 
momentum are then adjusted to reduce the eccentricity. 
The method introduced in [23|,[24| is capable of producing 
eccentricities of order of onl y 1 0~ 5 . For standard punc- 
ture data the method in [23ll24{ cannot be used directly. 
One problem is that standard moving puncture simula- 
tions start with the lapse and shift given in Eqs. (fT5|). 
Thus the coordinates used are not well adapted to quasi- 
equilibrium. Hence oscillations in the hole separation r 
are not due to real eccentricity alone but also due to the 
fact that the coordinates are still evolving as well. This 
problem is quite visible in Fig. [3] The broken line shows 
oscillations that have more than one frequency. Thus we 
cannot fit the curve very well with a straight line plus 
a single sine function as in (23l . l24j . To summarize, re- 
ducing the eccentricity for puncture data is harder than 
for the extended conformal thin sandwich data in [23| . 
Therefore, our results in both final eccentricities and re- 
duction of eccentricity per iteration are not as good as 

in mm. 

In the cases we have studied so far, we have found 
it unnecessary to adjust p r (away from the PN value in 
Eq. (pH]) ) in order to reach an eccentricity of order 0.001. 
It is much more important to choose an appropriate pt- 
However, we expect that adjusting p r will be necessary 
to reach even lower eccentricities. 



VI. DISCUSSION 

We have introduced the two new eccentricity measures 
e r and e22- Both are easy to compute since their calcu- 
lation does not involve any free parameters unlike earlier 
definitions akin to e u (t). To compute eccentricity mea- 
sures such as e w (t) one needs to specify a time interval 
during the inspiral phase over which we fit the orbital an- 



gular velocity to a polynomial of some low degree. This 
degree is essentially another free parameter and is usually 
chosen to be 4 or 5. Note, however, that all eccentricity 
definitions are ambiguous to a certain extend since the 
entire concept of eccentricity is only rigorously defined 
for periodic orbits. All eccentricity definitions for inspi- 
ral orbits depend on how we split a function of time (like 
the separation) into a smooth and an oscillating piece. 
Thus we do not claim that our definitions e r and e22 
are less ambiguous than earlier ones. For example our 
definitions depend on the period T which we compute 
from Kepler's law, but other choices are possible (e.g. 
an estimate of the actual orbital period of the numerical 
simulation at each time). Since our definitions e r and 
e22 do not require us to choose extra parameters such as 
fitting intervals they are easier to use. 

We also show that certain low resolution grid setups 
(which are relatively inexpensive) can be used to estimate 
the initial eccentricity using our measure e r . This gives 
us a relatively efficient way to measure the eccentricity 
of any initial data. 

Furthermore, we have explained why the coordinates 
of standard puncture data are not the same as in PN 
calculations in ADMTT or harmonic gauge. Thus even 
if one arrives at a highly accurate estimate of the mo- 
mentum parameters needed for low eccentricity orbits in 
ADMTT gauge, there is no easy way to incorporate this 
knowledge into simulations starting from standard punc- 
ture data. Indeed we find that using accurate parameters 
in ADMTT gauge [2l|, [22[ in standard puncture initial 
data can lead to relatively large eccentricities. 

We provide a simpler approach which starts from mo- 
mentum parameters using the relatively short expressions 
given by Kidder [4l|. After measuring the resulting ec- 
centricity e,. in a short inexpensive numerical simulation, 
these parameters can then be refined by changing p t by 
an amount of order Ap t (see Eq. ([23|) 1. In this way we 
can always arrive at eccentricities that are low enough for 
the purposes of the NRAR collaboration 19]. In order 
to achieve this we usually do not need more than three 
numerical runs. 



Acknowledgments 

It is a pleasure to thank Doreen Miiller for use- 
ful discussions about the approach in [22j . This work 
was supported by NSF grant PHY-0855315. Compu- 
tational resources were provided by the Ranger cluster 
at the Texas Advanced Computing Center (allocation 
TG-PHY090095) and the Kraken cluster (allocation TG- 
PHY100051) at the National Institute for Computational 
Sciences. 



[1] B. Abbott et al. (The LIGO Scientific Collaboration) [2] LIGO, http://www.ligo.caltech.edu. 
(2007), arXiv:0711.3041 [gr-qc], 0711.3041. 



8 



[3] F. Acernese et al., J. Opt. A: Pure Appl. Opt. 10, 064009 
(2008). 

[4] VIRGO, http://www.virgo.infn.it. 

[5] GEO600, http://www.geo600.org. 

[6] B. Schutz, Class. Quantum Grav. 16, A131 (1999). 

[7] T. Nakamura, K. Oohara, and Y. Kojima, Prog. Theor. 

Phys. Suppl. 90, 1 (1987). 
[8] M. Shibata and T. Nakamura, Phys. Rev. D 52, 5428 

(1995). 

[9] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D59, 
024007 (1998), gr-qc/9810065. 
[10] W. Tichy, Phys. Rev. D74, 084005 (2006), gr- 
qc/0609087. 

[11] W. Tichy, Class. Quant. Grav. 26, 175018 (2009), 
0908.0620. 

[12] W. Tichy, Phys. Rev. D80, 104034 (2009), 0911.0973. 
[13] M. Campanelli, C. O. Lousto, P. Marronetti, and 

Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr- 

qc/0511048. 

[14] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr- 
qc/0511103. 

[15] S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 

(1997), gr-qc/9703066. 
[16] M. Ansorg, B. Briigmann, and W. Tichy, Phys. Rev. 

D70, 064011 (2004), gr-qc/0404056. 
[17] P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963). 
[18] P. C. Peters, Phys. Rev. 136, B1224 (1964). 
[19] NRAR, https://www.nmja- 

pro j ect . org/doku . php ? id = nrar : home . 
[20] J. G. Baker, J. R. van Meter, S. T. McWilliams, J. Cen- 
trella, and B. J. Kelly, Phys. Rev. Lett. 99, 181101 

(2007), gr-qc/0612024. 
[21] S. Husa, M. Hannam, J. A. Gonzalez, U. Sperhake, 

and B. Bruegmann, Phys. Rev. D77, 044037 (2008), 

0706.0904. 

[22] B. Walther, B. Bruegmann, and D. Mueller, Phys. Rev. 

D79, 124040 (2009), 0901.0993. 
[23] H. P. Pfeiffer et al., Class. Quant. Grav. 24, S59 (2007), 

gr-qc/0702106. 

[24] M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mrou, H. P. 
Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, 
Phys. Rev. D76, 124038 (2007), 0710.0158. 

[25] A. H. Mroue, H. P. Pfeiffer, L. E. Kidder, and S. A. 
Teukolsky (2010), 1004.4697. 

[26] A. Buonanno, G. B. Cook, and F. Pretorius, Phys. Rev. 



D75, 124018 (2007), gr-qc/0610122. 
[27] B. Briigmann, W. Tichy, and N. Jansen, Phys. Rev. Lett. 

92, 211101 (2004), gr-qc/0312112. 
[28] B. Briigmann, J. Gonzalez, M. Hannam, S. Husa, 

U. Sperhake, and W. Tichy, Phys. Rev. D77, 024027 

(2008), gr-qc/0610128. 
[29] P. Marronetti, W. Tichy, B. Briigmann, J. Gonzalez, and 

U. Sperhake, Phys. Rev. D77, 064010 (2008), 0709.2160. 
[30] M. J. Berger and J. Oligcr, Journal of Computational 

Physics 53, 484 (1984). 
[31] T. W. Baumgarte, Phys. Rev. D 62, 024018 (2000), gr- 

qc/0004050. 

[32] W. Tichy, B. Briigmann, and P. Laguna, Phys. Rev. 

D68, 064008 (2003), gr-qc/0306020. 
[33] W. Tichy and B. Briigmann, Phys. Rev. D 69, 024006 

(2004), gr-qc/0307027. 
[34] P. Marronetti, W. Tichy, B. Briigmann, J. Gonzalez, 

M. Hannam, S. Husa, and U. Sperhake, Class. Quant. 

Grav. 24, S43 (2007), gr-qc/0701123. 
[35] B. Briigmann, J. A. Gonzalez, M. Hannam, S. Husa, 

U. Sperhake, and W. Tichy, Phys. Rev. D77, 024027 

(2008), gr-qc/0610128. 
[36] G. Schafer, Annals of Physics 161, 81 (1985). 
[37] W. Tichy, B. Briigmann, M. Campanelli, and P. Diener, 

Phys. Rev. D67, 064008 (2003), gr-qc/0207011. 
[38] N. Yunes and W. Tichy, Phys. Rev. D74, 064013 (2006), 

gr-qc/0601046, gr-qc/0601046. 
[39] N. Yunes, W. Tichy, B. J. Owen, and B. Briigmann, 

Phys. Rev. D74, 104011 (2006), gr-qc/0503011. 
[40] N. K. Johnson-McDaniel, N. Yunes, W. Tichy, and B. J. 

Owen, Phys. Rev. D80, 124039 (2009), 0907.0891. 
[41] L. E. Kidder, Phys. Rev. D52, 821 (1995), gr- 

qc/9506022. 

[42] W. Tichy and P. Marronetti, Phys. Rev. D76, 061502 

(2007) , gr-qc/0703075. 

[43] W. Tichy and P. Marronetti, Phys. Rev. D78, 081501 

(2008) , 0807.2985. 

[44] W. Tichy, B. Briigmann, and P. Laguna, Phys. Rev. 

D68, 064008 (2003), gr-qc/0306020. 
[45] J. W. York, Phys. Rev. Lett. 82, 1350 (1999). 
[46] The standard PN jargon for a term of order (v/c) n in the 

equations of motion is ^PN term. Notice that according 

to the virial theorem v 2 ~ M/r, so that ^PN terms are 

also 0(M/r)%. 



